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Abstract 

Using molecular dynamics simulations and integral equations (Rogers- Young, Percus-Yevick and 
hypernetted chain closures) we investigate the thermodynamic of particles interacting with contin- 
uous core-softened intermolecular potential. Dynamic properties are also analyzed by the simula- 
tions. We show that, for a chosen shape of the potential, the density, at constant pressure, has a 
maximum for a certain temperature. The line of temperatures of maximum density (TMD) was 
determined in the pressure-temperature phase diagram. Similarly the diffusion constant at a con- 
stant temperature, D, has a maximum at a density Pmax aud a minimum at a density Pmin < Pmax- 
In the pressure-temperature phase-diagram the line of extrema in diffusivity is outside of TMD 
line. Although in this interparticle potential lacks directionality, this is the same behavior observed 
in SPC/E water. 

PACS numbers: 64.70.Pf, 82.70.Dd, 83.10.Rs, 61.20.Ja 
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I. INTRODUCTION 



Water is an anomalous substance in many respects. Most liquids contract upon cooling. 
This is not the case of water, a liquid where the specific volume at ambient pressure starts 
to increase when cooled below T = A°C Besides, in a certain range of pressures, also 

□ n 

exhibits an anomalous increase of compressibility and specific heat upon cooling Far 
less known are its dynamics anomalies: while for most materials diffusivity decreases with 
increasing pressure, liquid water has an opposite behavior in a large region of the phase 
diagram [3, , Q, S Is , [lO , U , 12 , Q . The increase of diffusivity of water as the pressure is 
increased is related to the competition between the local ordered tetrahedral structure of the 
first neighbors and the distortions of the structure of the first and second neighbors. In the 
region of the phase diagram where this ordered structure is dominant, increasing pressure 
implies breaking first neighbors hydrogen bonds what allow for interstitial second neighbors 
to be in a closer approach. The interactions are thus weakened and therefore, although 
the system is more dense it has a larger mobility. In this sense, a good model for water 
and tetrahedral liquids should not only exhibit thermodynamic but also dynamic anomalies. 
In SPC/E water, the region of the pressure-temperature (p -T) phase-diagram where the 
density anomaly appears is contained within the region of the p -T phase-diagram where 
anomalies in the diffusivity are present 0, lo| . 

For explaining the thermodynamic anomalies, it was proposed that these anomalies are 
related to a second critical point between two liquid phases, a low density liquid (LDL) 
and a high density liquid (HDL) [l^ located at the supercooled region beyond the fine of 
homogeneous nucleation and thus it cannot be experimentally measured. 

Water, however, is not an isolated case. There are also other examples of tetrahedrally 
bonded molecular liquids such as phosphorus Q and amorphous silica |l7 1 that also are 
good candidates for having two liquid phases. Moreover, other materials such as liquid metals 
jisl and graphite also exhibit thermodynamic anomalies. Unfortunately a coherent and 
general interpretation of the low density liquid and high density hquid phases is still missing. 

What type of potential would be appropriated for describing the tetrahedrally bonded 
molecular liquids? Directional interactions are certainly an important ingredient in obtaining 
a quantitative predictions for network-forming liquids like water. However, the models that 
are obtained from that approach are too complicated, being impossible to go beyond mean 
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field analysis. Isotropic models became the simplest framework to understand the physics of 
the liquid-liquid phase transition and liquid state anomalies. From the desire of constructing 
a simple two-body isotropic potential capable of describing the complicated behavior present 
in water-like molecules, a number of models in which single component systems of particles 
interact via core-softened (CS) potentials have been proposed. They possess a repulsive 



core that exhibits a region of softening where t 



can be a shoulder or a ramp 
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y. This region 



311 132, 1331 |3j,|35|,|3f 



In the first case, the potential consists of a hard core, a square repulsive shoulder and, in 



2^, 22,Q,HQ,QQ• 



In two 



some cases, an attractive square well 
dimensions, such potentials have density and diffusion anomalies and in some cases a second 



critical point 
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In three dimensions, these potentials do not have dynamic 
and thermodynamic anomalies but possess a 
point, access! 
equation [2^ 



ible by simulations in the region predicted by the hypernetted chain integral 

In the second case, the interaction potential has two competing equilibrium distances, 
defined by a repulsive ramp 0, Q, Q. By including a global term for attractions, this 
model displays a liquid phase with a first-order line of liquid-gas transition ending in a 



critical point 



33l l36| brings this second critical point in to an accessible region of higher 
temperature, and also displays a normal gas-liquid critical point. 

Notwithstanding the progresses made by the models described above, a potential in which 
both the potential and the force are continuous functions and that exhibits all the thermo- 
dynamic and dynamic anomalies present in water is still missing. In this paper, we check 
if particles interacting with a core-softened potential similar to the one proposed by Cho et 



al. 



37 



38, 



39| exhibit thermodynamic and dynamic anomalies similar to the ones present 



in water. Since the potential can have a variety of shapes, depending on its parameters, we 
study a soft ramp ( with continuous force) with two scale distances. This tvpe of potential 
gives a distribution function similar to the one expected for SPC/E water 40|. We check if 
the region in the pressure-temperature phase-diagram of thermodynamic anomalies is inside 
the region of dynamic anomalies as in SPC/E water Iic(. 

The reminder of this paper goes as follows. In sec. II the model is introduced; in sec. 
Ill the phase-diagram is obtained within the Rogers- Young, Percus-Yevick and hypernetted 
chain integral equations. Results for the phase-diagram and for the diffusion constant ob- 
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tained from molecular dynamics simulations are shown in sec. IV. Conclusions about the 
relation between the locus of the density anomaly and the diffusion anomaly are presented 
in sec. V. 



II. THE MODEL 



We consider a set of molecules of diameter a interacting through a potential that consists 
of a combination of a Lennard- Jones potential of well depth e plus a Gaussian well centered 
on radius r = tq with depth a and width c, 

2" 



f/(r) = 4e 



12 



(T\6 



r — To 
a 



(1) 



This potential can represent a whole fami 
tions, from a deep double wells potential 0, 



of two length scales intermolecular interac- 
^ to a repulsive shoulder j^^, depending on 



the choice of the values of a, ro and c. Specific choices of these parameters leads to double 



well potentials similar to the one studied by Cho et al. [37[. The attractive double well 
brings both the liquid-gas phase-transition and the anomalies to higher temperatures into 
the unstable region of the p -T phase diagram 3^. 

In order to circumvent this difficulty, here we investigate the thermodynamic and dynamic 
behavior of particles interacting via a potential with a very small attractive region. We use 
Eq. (0) with a = 5, r^/a = 0.7 and c = 1. This potential has two length scales within a 
repulsive ramp followed by a very small attractive well ( Fig. (Q) ). 

In order to have an overview of the behavior of particles interacting with this potential, 
we use integral equations to estimate the thermodynamic properties in the phase diagram. 



III. INTEGRAL EQUATIONS 

One of the most successful theories for describing the structure of simple fluids are the 
integral equations Among them, certainly the most famous is the Ornstein-Zernike 

(OZ) equation j^] that, for pure isotropic fluids with density p, gives an exact relation 
between the direct correlation function, c(r), and the total correlation function, h(r), and it 
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FIG. 1: Interaction potential eq. with parameters a = 5, ro/cr = 0.7 and c = 1, in reduced 
units. The inset shows a zoom in the very small attractive part of the potential 

is given by 

7(r) = h(r) — c(r) ~ P J h(r)c(\r — r \)dr, (2) 

where h{r) = g{r) — l and where g{r) is the pair distribution function. g{r) is proportional 
to the probability to find a particle at a distance r when another particle is placed at the 
origin. 

The Fourier transform of Eq. Q is given by 

^^"'-l- pC{k) ' 

where T{k) and C{k) are the Fourier Transforms of 7(r) and c(r), and the definition for the 
direct correlation function 

c(r) = h{r) - In {g{r) exp [(3U{r)]} + B{r), (4) 

was used. Here (3 = l/ksT and B{r) is the sum of all bridge diagrams for the interparticle 
potential. Eq. Q together with Eq. (jH) can be solved for a given interparticle potential. 
For obtaining Bir) many approximations {closure relations) have been proposed 
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48l |49| along the years. Unfortunately, these approximations have the following 
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FIG. 2: Pressure-temperature phase-diagram obtained by Rogers- Young integral equations. 
From bottom to top, the isochores p = 0.100; 0.103; 0.105; 0.107; 0.110; 0.113; 0.115; 0.117; 0.120; 
0.123; 0.125; 0.127; 0.130; 0.132; 0.134; 0.136; 0.138; 0.140; 0.142 and 0.144 are shown. The sohd line 
illustrate the TMD. 

thermodynamic inconsistence: the pressure calculated via the fluctuations route, 





differs from the pressure calculated via the virial route, 




(6) 



Two of these closures have been widely used: the Percus-Yevick (PY) 




where 



B{r) = ln[l + 7(r)] - 7(r) , 



(7) 



and the hypernetted chain (HNC)|4 



that sets 



B{r) = . 



(8) 
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FIG. 3: Phase-diagram in Rogers- Young integral equations. The low isochores p = 
0.01; 0.02; 0.03; 0.04; 0.05; 0.06; 0.07 and 0.08, from bottom to top, are shown. The inset shows 
the isochores for p* < 0.03 converging to a point at supercooled region - the liquid-gas phase 
transition. 

While HNC is appropriated for large interparticle distances, PY is more adequate for 
small ones. In order to avoid the inconsistencies present in the original integral equations, 
Rogers and Young js^ proposed a mix of the HNC and PY closures of the form 



with the mixing function /(r) = 1 — exp[— ar]. Note that at r = Eq. Q reduces to 
the PY approximation and for r cxd, Eq. © tends to the HNC approximation. The 
Rogers- Young (RY) approximation puts together PY and HNC closures with an adjustable 
parameter a. This parameter is determined by imposing that the pressure calculated using 
Eq. © gives the same result as using Eq. (jH)) [global consistency criterion). This method 
have the inconvenience of the integral in p' . Instead of calculating a by imposing that the 
pressure should be the same when calculated using Eq. (0) and Eq. one can obtain 
a by checking the consistency between the compressibilities Xfiuc and Xvin calculated by 



c(r) = exp[-/5[/(r)] 1 + 



exp[7(r)/(r)] - 1 
f{r) 



7(r) - 1 



(9) 
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derivation of Eqs. (0) and Q respectively j5l|, namely 



P ' ' Jo 

and 



3 C"" 
XfL=^-^^P r\{r)dr (10) 



27T ^ r .dU{r)dg{p,r) ^ 



52] or more adjustable parameters 



53, 



Others closures were proposed where one j51| 

y, 

are needed in order to guarantee the consistency. In this work, we use the RY 
approximation due its success in describing the structure of the systems whose particles 
interact by a purely repulsive pair potentials 351 l5Cll. l56l l57|. 



A numerical iterative solution of the system formed by Eq. (jH)) and Eq. Q was performed 
using a fine grid with M = 4096 points and a step size Ax = 0.0075, from x = r/a = 0.0075 
until X = MAx. The tolerance for thermodynamic consistency was 1 — Xfiuc./Xvir. < 10^^. 
For the PY and HNC closures, the same M and Ax was used, in the same range. Pressure, 
temperature and density are shown in dimensionless units: 

T* ^ M (12) 



P* = per' (13) 

P* = t-^ (14) 

The main features of the phase diagram obtained by RY closure are illustrated in Fig. 
This p -T phase-diagram shows that the isochores with 0.120 < p* < 0.140 have minimum 
which means that {dP/dT)p = 0. From this follows {dp/dT)p = 0, which implies a density 
anomaly. The line of minima for the different isochores forms the Temperatures of Maximum 
Density (TMD) shown in Fig. (0) by a solid line. 

The presence of a possible critical point between two liquid phases may 
be suggested by the crossing of the analytic continuation of isochores p* = 
0.134; 0.136; 0.138; 0.140; 0.142; 0.144, in the region below T* < 0.05. In this region the 
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integral equations numerical solutions do not converge and the thermodynamic equilibrium 
is not achieved for the MD simulations. 

Fig. © shows that the very low isochores at p* < 0.03 are converging to a point in the 
supercooled region, indicating a liquid-gas critical point, as can be seen from the inset. 

When analyzing the model Eq. (Q) with the PY approximation, density anomaly was 
found between 0.13 < p* < 0.3 in a region of temperatures from 0.4 < T* < 0.86, as can be 
seen from Fig. (@^). No density anomaly was found when the model Eq. (0) was analyzed 
with HNC (Fig. (gb)). 




0.3 0.6 0.9 1.2 0.3 0.6 0.9 1.2 

T* T* 

FIG. 4: p -T phase-diagram obtained by PY (Fig.(@t)) and HNC (Fig. (jit)) integral equations. 
From bottom to top, the twenty five isochores illustrated are p = 0.10; 0.11; 0.12; 0.33; 0.34. in 
both figures. The solid line in Fig. (a) illustrates the TMD line. 

IV. THE MOLECULAR DYNAMICS 

We also performed molecular dynamics simulations in the canonical ensemble using 500 
particles in a cubic box with periodic boundary conditions, interacting with the intermolecu- 
lar potential described above. The chosen parameters were a = 5.0, r^/a = 0.7 and c = 1.0. 
The cutoff radius was set to 3.5 length units. Using reduced units defined as T* and p*, 
a broad range of temperatures (0.10 < T* < 0.45) and densities (0.05 < p* < 1.00) was 
chosen, in order to explore the phase diagram. Thermodynamic and dynamic properties 
were calculated over 2 500 000 steps long simulations, previously equilibrated over 500 000 
steps. In the lower temperature systems, additional simulations were carried out with equi- 
libration over 2 000 000 steps, followed by 6 000 000 simulation run. The time step was 
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FIG. 5: Pressure-temperature phase-diagram obtained by molecular dynamic simulation. From 
bottom to top, the same isochores illustrated in RY phase-diagram, p = 0.100; 0.103; 0.105; 0.107 
0.110; 0.113; 0.115; 0.117; 0.120; 0.123; 0.125; 0.127; 0.130; 0.132; 0.134; 0.136; 0.138; 0.140; 0.142 and 
0.144 are shown. The solid line illustrates the TMD and the dashed line shows the boundary of 
the diffusivity extrema. 

0.001 in reduced units. The thermodynamical stability of the system was checked analyzing 
the dependence of pressure on density and also by visual analysis of the final structure, 
searching for cavitation. 

Figure (0) shows the p -T phase-diagram obtained by molecular dynamics. The isochores 
have minima that define the temperature of maximum density. The TMD line encloses 
the region of density (and entropy) anomaly. The comparison between the RY and MD 
results shows that the TMD line starts at lower densities in the MD simulations than that 
in RY integral equations. Above p* = 0.144 both theories agree that no density anomaly 
happens. The RY pressures for each isochore are slightly underestimated when compared 
with simulations, but the overall agreement between the predictions of this closure and the 
simulations results is very good. On the other side, the PY approximation predicts density 
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FIG. 6: Diffusion coefficient as a function of density for some studied temperatures. Tlie units 
are defined in tlie text. 

anomaly, but in a region completely different than that of MD. The HNC closure do not 
shows TMD line, as we have discussed before. 

The MD simulations also indicates the possibility of a liquid-liquid critical point by the 
crossing of the analytic continuation of isochores p* = 0.134; 0.136; 0.138; 0.140; 0.142; 0.144 
- the same behavior was seen by the RY closure, and missing by PY and HNC. 

We also study the mobility associated with the potential described in Eq. (jT)). The 
diffusion is calculated using the the mean-square displacement averaged over different initial 
times, 

(Ar(tn = ([r(to + t)-r(toF). (15) 
Then the diffusion coefficient is obtained from the relation 



D = lim(Ar(t)2)/6t . 

t— »oo 



Figure (jH} shows the behavior of the translational diffusion coefficient, 

D(m/e)i/2 



D* 



a 



(16) 



(17) 
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as a function of p* . At low temperatures, the behavior is similar to the behavior found [lO| 
in SPC/E supercooled water. The diffusivity increases as the density is lowered, reaches a 
maximum at p_Dmax (and Pomax) and decreases until it reaches a minimum at porain (and 

-Pomin) • 

The region in the p -T plane where there is an anomalous behavior in the diffusion is 
bounded by (Tomin, -Pomin) and (Tomax, -Pomax) and their location is shown in Fig. (0). The 
region of diffusion anomalies (Tomax, PVwx) and (Tomin, pDmin) lies outside the region of 
density anomalies like in SPC/E water [10|. 



V. CONCLUSIONS 



We have studied the thermodynamic properties of fluids interacting via a three dimen- 
sional continuous core-softened potential with a continuous force, using several integral equa- 
tions closures, RY, FY and HNC, as well as molecular dynamics simulations. The continuity 
of the force is similar that one expect for realistic systems. We studied the density anomaly 
and anomalies in the translational diffusion. Both RY integral equations and molecular 
dynamics results show that the density can behave anomalously at a certain range of pres- 
sures and temperatures. The agreement between these two theories is very good, confirming 
the RY integral equations as a powerful tool for investigations of interatomic repulsive pair 
potentials. The FY approximation emphasizes the short-ranged interactions and indeed 
predicts density anomaly, but the width of the anomaly region is strongly overestimated if 
compared with MD simulations. No density anomaly was found employing the HNC closure, 
because this approach is better suited for systems with long-ranged interactions. 

Both MD and RY theories suggests the possibility of a second critical point, between two 
liquid phases, by the crossing of the analytic continuation of the isochores where 0.134 < 
p* < 0.144 for T* < 0.05. However the actual calculations or simulations in this region 
was not possible, either by failure in the integral equations convergence and because the 
equilibrium was not reached. 

The translational diffusion shows a maximum and a minimum in the pressure-temperature 
phase-diagram. The region in the p -T plane of density anomaly is located inside the region 
of the anomalous diffusion. 

The studied continuous core-softened potential, despite of not having long-ranged or H- 
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bond-like directional interactions, exhibit thermodynamic and dynamic anomalies similar to 
the ones observed in SPC/E water J^]. 
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